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ABSTRACT 

We perform population synthesis studies of different types of neutron stars (thermally 
emitting isolated neutron stars, normal radio pulsars, magnetars) taking into account 
the magnetic field decay and using results from the most recent advances in neutron 
star cooling theory. For the first time, we confront our results with observations using 
simultaneously the Log N - Log S distribution for nearby isolated neutron stars, the 
Log N - Log L distribution for magnetars, and the distribution of radio pulsars in the 
P - P diagram. For this purpose, we fix a baseline neutron star model (all microphysics 
input), and other relevant parameters to standard values (velocity distribution, mass 
spectrum, birth rates ...), allowing to vary the initial magnetic field strength. We find 
that our theoretical model is consistent with all sets of data if the initial magnetic 
field distribution function follows a log-normal law with < log(Bo/[G]) >^ 13.25 
and oi og B ~ 0.6. The typical scenario includes about 10% of neutron stars born 
as magnetars, significant magnetic field decay during the first million years of a NS 
life (only about a factor of 2 for low field neutron stars but more than an order 
of magnitude for magnetars), and a mass distribution function dominated by low 
mass objects. This model explains satisfactorily all known populations. Evolutionary 
links between different subclasses may exist, although robust conclusions are not yet 
possible. 
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1 INTRODUCTION 

At the present moment our knowledge about neutron star 
(NS) evolution is an intriguing puzzle. We know many ob- 
servational manifestations of young isolated NSs: radio pul- 
sars (PSRs); central compact objects in supernova rem- 
nants (CCOs in SNRs); rotating radio transients (RRATs); 
radio-quiet thermally emitting isolated NSs, also known 
as X-ray dim isolated NSs (XDINS) or the Magnificent 
Seven (M7), and the observational manifestations of mag- 
netars -soft gamma-ray repeaters (SGRs) and anomalous 
X-ray pulsars (AXPs). Reasons for this apparent diver- 
sity as well as possible links between different classes 
are not entirely clear (see e.g. Manchester et all 
Woods fc Thompson! ll200r3): iKaspil il2007h ; lHaberl 



2005) 



2007) 



Zand (|2007h ; iRea fe McLaughlin! (|2008h for recent reviews 
about the different subclasses). 

In the past few years we have learned that some 
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NSs can show different types of activity, transiting from 
class to class. For example, PSR J1846-0258, known for 
some time as a normal pulsar ; demonstrated outbursts 
typical for AXPs o r/and SGRs jKumar fc Safi-Harbl liooj ; 
iGavriil et al. 2008). In addition, the total energy release 
by the object became larger than the rotational energy 
losses, w hich according to the origin al classification dis- 
cussed in iThompson fc Duncan! (|l996l ) should place it in 
the AXP list. Thus, for the first time, we have an exam- 
ple of transformation of a PSR into a magnetar. Several 
of the SGRs do not show any bursting activity for many 
years, and if we had not enough information about their vi- 
olent past, we would have classified them as AXPs. Some 
of R RATs were shown to emit normal radio p ulsar emis- 
sion (|Deneva et aill2008l : iMcLaughlin et~aH2009l h The tran- 
sient AXP XTE J181 0-197 and AXP IE 1547.0- 5408 demon- 
strated radio pulses (ICamilo et alj|2006l 120071 ). One of the 
RRATs J1 819-1458 shows therm al properties much similar 
to the M7 l|Revnolds et al.ll2006l ). Hence, divisions between 
some subpopulations of young isolated NSs (or at least some 
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of their representatives) can be illusive. On the other hand, 
young, low-field CCOs (jHalpern et alj|2007h . normal PSRs 
(B ~ 10 12 G, i.e. Crab and Vela-like), and SGRs clearly 
represent NSs born with different properties. 

Our brief observational record of NSs (~ 40 years at 
most), low statistics in many cases, and selection effects do 
not allow to draw a coherent picture of NS "sociology" just 
from observations. From the theoretical point of view, our 
understanding of the SN explosion mechanism is not pre- 
cise enough to provide a solid model of initial parameters of 
NSs and evolutionary models (thermal, magnetic field, and 
spin evolution) are related to extremely complicated physi- 
cal problems (superfluidity and superconductivity in dense 
matter, electrodynamics in superstrong magnetic fields, etc.) 
which usually leads to inconclusive results. 

In our opinion, it is necessary not only to compare ob- 
servations with models to verify individual objects, but also 
to confront theoretical calculations with observational data 
via population synthesis techniques taking into account as 
many classes of NSs as possible. Joint constraints by means 
of simultaneous comparison of theoretical models with dif- 
ferent subpopulations should be derived to form a population 
mosaic. Numerous degrees of freedom in modern evolution- 
ary models must be compensated by different observational 
tests. Several important s tudies in this area appeared in re - 
cent years, see for examp le lFaucher-Giguere fc Kaspll (2006) ; 
iKeane fc Kramerl (|2008l ). and references therein. In this pa- 
per we continue to follow this lane but with an important 
difference: we attempt to constrain our model and check its 
consistency using at the same time different populations. 
This is, to our knowledge, the first time that a multilat- 
eral approach is employed. Previous works that focused in 
a specific NS observational class, ignored that the param- 
eters obtained may be in contradiction with the properties 
of a different class of objects. For examples, models without 
magnetic field decay are clearly in contradiction with the 
existence of magnetars. We try to give a step forward in the 
direction towards a NS unified model by using at the same 
time different populations. 

Among the different physical ingredients needed to 
properly model the thermal evolution of NSs, we empha- 
size that heat transport in the NS crust plays a crucial 
role during the first millio n years of its thermal evolution 
jYakovlev fc Pethickl liooi) . which is typically the period 
during which NSs are detectable with current X-ray instru- 
ments. Multi-wavelength observations in the soft X-ray, UV, 
and optical bands of the thermal emission from a NS's sur- 
face now provide a real opportunit y to probe the internal 
physics of NSs fsee lPage et"ai1l2006l for a general overview) . 

Remarkably, all isolated nearby compact X-ray sources 
that have been detected also in the optical band (RXJ 
185635-3754, RX J0720.4-3125, RX J1308.6+2127, and RX 
J1605. 3+3249) have a significant optical exc ess relative t o 
the extrapolated X-ray blackbody emission l|Haber]||2007h . 
An inhomogeneous surface temperature distribution can 
accommodate this optical excess and can arise naturally 
if heat conduction in the NS crust is a nisotropic due 
to the presence of a large magnetic field ijGeppert et al.l 
|2004 120061 ; IPerez-Azorfn et al.ll2006l ). Alternative models 
are based on anisotropic radi ation from magnetized atmo- 
spheres, as in iHo et all J2007h . or condensed surfaces as in 
IPerez-Azorfn et al.l ( 20051 ). None of the explanations (only 



temperature anisotropy vs. physical origin) are entirely sat- 
isfactory and probably a combination of both effects is 
needed. Heat conduction can also influence other observ- 
able aspects of accreting NSs in low mass X-ray binaries, 
including their quiesc ent luminosity and the superburst re- 
currence time-scales l|Brown et al.| [l998l ). For this reason, 
one of the goals of this paper is to revisit former population 
synthesis studies using new NS evolution models that in- 
clude anisotrop ic heat transport in N Ss crusts and magne tic 
field evolution (lAguilera et al.ll2008al lbl; iPons et al.ll2009l ). 

In this paper we want to study if the patchy view of 
young NSs subpopulations can be explained by a unique 
set of smooth distributions of the most important parame- 
ters, among which the the magnetic field distribution plays 
the main roleQ Our main idea here is to use several dif- 
ferent tests to confront theoretical predictions and observa- 
tions. For this purpose, the paper is organized as follows. In 
the next section we briefly describe the model of magneto- 
thermal evolution of NSs that we use. With this model, one 
has the cooling behavior and magnetic field (and therefore 
period) evolution of NS. In Sec. 3, we describe the popula- 
tion synthesis technique used for Log N - Log S calculations 
of close-by cooling NSs which thermal emission has been 
detected and which temperature has been estimated. With 
this, we can partially constrain the initial magnetic field 
distribution. Then, in Sec. 4, we discuss the Log N - Log L 
distribution of the population of magnetars in the Galaxy. 
We show that the model constrained in the previous section 
is also consistent with the flux distributions of the extrapo- 
lated magnetar population. To finish the presentation of our 
results, in Sec. 5 we use the population of radio pulsars in the 
P — P diagram to put additional constraints on the proper- 
ties of NSs. Explicitly, the current period and magnetic field 
distributions of rotation-powered pulsars constrain the NS 
initial magnetic field distribution and break the degeneracy 
in the parameter space obtained in previous section. Sec. 6 
is devoted to the final remarks and to discuss uncertainties 
of the models we use and future prospects. 



2 MAGNETIC FIELD DECAY AND COOLING 
MODEL 

Very often thermal and magneto-rotational evolution of NSs 
are treated separately. However, in the case of young (< 1 
Myr) NSs with magnetic fields > 10 13 G this is incorrect, be- 
cause temperature affects the electrical resistivity, and there- 
fore the magnetic field evolution, while field decay provides 
an additional energy source that modifies the temperature 
of the star. Although for the average radio pulsar population 
(old and rel atively low field) this effect is pr obably not very 
important |Faucher-Giguere fc Kaspll |2006i ), there is some 
observational evidence of the interplay of the magnetic field 
and temperature during early stages of NS evolution. As 



Note for the arXiv version: After the final (second revision) ver- 
sion of this paper has been submitted, an important e-print by 
Kaplan and van Kerkwijk appeared (arXiv: 0909.5218). In this 
article the authors give new observational arguments in favor of 
the large role of the field decay in the evolution of the Magnifi- 
cent Seven. Many of our conclusions coincide with those given by 
Kaplan and van Kerkwijk. 
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discussed in iPons et all (|2007l ), there is a strong correlation 
between the inferred magnetic field and the surface temper- 
ature in a wide range of magnetic fields: from magnetars 
(B > 10 14 G), through radio-quiet isolated NSs (B ~ 10 13 
G) down to some ordinary PSRs (B < 10 13 G). The main 
conclusion was that, rather independently from the stellar 
structure and the matter composition, the correlation can 
be explained by heating from dissipation of currents in the 
crust on a timescale of ~ 10 6 yrs. 

This observed correlation has been confirmed later by 
more detailed 2D cooling simulations combining the insu- 
lating effect of strong non radial fields with the additional 
source of heating due the Ohmic dissipation of the ma g- 
netic field in th e crust al region l|Aguilera et all l2008al lbl; 
IPons et al.ll2009t ). It was shown that, during the neutrino 
cooling era and the early stages of the photon cooling era, 
the feedback between Joule heating and magnetic diffusion is 
strong, resulting in a faster dissipation of the stronger fields. 
As a consequence, all neutron stars born with fields over a 
critical value (> 5 x 10 13 G) reach similar field strengths 
(~ 2 — 3 x 10 13 G) at late times. Irrespective of the initial 
magnetic field strength, the temperature becomes so low af- 
ter a few million years that the magnetic diffusion timescale 
becomes longer than the typical ages of PSRs, thus appar- 
ently resulting in no dissipation of the magnetic field in old 
NS. Another interesting result was that the effective temper- 
ature of models with strong internal toroidal components 
is systematically higher than that of models with purely 
poloidal fields, due to the additional energy reservoir stored 
in the toroidal field that is gradually released as the field 
dissipates. 

In this paper, we have employed cooling curves obtained 
fro m an upda t ed ver sion of the 2D cooling code described 
m IPons et al l (|2009l ), that also includes the effect of su- 
perflu id heat conduction in the inner crust (|Aguilera et al.l 
I2009T ). We have used a Skyrme-type equation of state (EoS) 
at zero temperature describing both, the NS crust and the 
liquid core, based on the effective nuclear interaction SLy 
(|Douchin fe Haensell [200ll V It is a purely hadronic, rela- 
tively stiff EOS that gives typical proper radii in the range 
~ 11.3 — 11.8 km (the radius observed at infinity would 
be 10-30 % larger depending on the mass), while the crust 
thickness varies from 0.6 to 1.2 km also depending on the 
mass of the NS. 

We refer to section 4 in lAguilera et ail (|2008ah for 
more details about the cooling models (neutrino emissivi- 
ties, equation of state, thermal conductivities, etc.). Details 
about the magne tic field initial geometry can be found in 
IPons et al.l l|2009l ). Our baseline initial model consists of a 
crustal-confined magnetic field with a poloidal component, 
parameterized by the value of the radial component at the 
magnetic pole (B) combined with a toroidal component with 
a maximum v a lue of twice B (see Eqs. (11) and (13) of 
lAguilera et al.l ((20083)) . We find that, although the ampli- 
tudes of both fields are of the same order of magnitude, the 
contribution of the toroidal field to the total magnetic en- 
ergy is < 10%, because this field is non vanishing only in a 
finite region of the star. This model is in agreement with the 
results obtaine d in recent studies of magnetic equilibr ium 
configurations (|Ciom et al.ll2009l ; lLander fc Jonesil2009h . 

In the left panel of Fig. [TJ we show a sample of cool- 
ing curves (effective temperature versus true age) for an 



M = 1.25M© NS but varying the initial strength of the 
magnetic field. For B > 10 13 G the presence of strong mag- 
netic fields has a visible effect from the very beginning of 
the evolution (results for a M = IAMq are very similar to 
this case). The effective temperature of a young, t — 10 3 
yr magnetar with B > 10 15 G is a few times higher than 
that of a NS with a standard B — 10 13 G, and it is kept 
above 10 6 K for a much longer time. The effect is quali- 
tatively similar, although somewhat smaller, for high mass 
stars, and of course it also depends on details about the mi- 
crophysics input. The influence of the superfluid neutron gap 
in the core is particularly relev ant, in this work we use the 
results from lBaldo et all (|l998l h We also use an updated 1 5o 
neutron superfluid gap in the crust obtained from Quantum 
Monte Carlo simulations l|Gezerlis fc Ca rlson 2008|) . This re- 
sults in small diff erences when comparing to the results in 
IPons et al.ll|200g|) . In non-magnetized NSs, a smaller gap re- 
sults in higher temperatures at early times (suppression of 
neutrino emissivity), but varying the gap does not change 
significantly the temperature of magnetars. The purpose of 
this paper is to study the observational implications of dif- 
ferent initial magnetic fields, being all the rest equal, so 
hereafter we will fix the underlying physical model and we 
will only vary the normalization of the field (not its geom- 
etry) and the mass of the NS. The variability of the cool- 
ing curves on the mass for a fixed magnetic field is shown 
in the right panel of Fig. [T] We have chosen B — 10 14 G 
and varied the mass of the NS in a wide range, namely 
M = 1.10, 1, 25, 1.32, 1.40, 1.48, 1.60, 1.70, and 1.76 M©. Ex- 
cept for the two lowest masses, the rest of cooling curves are 
difficult to be distinguished. This is clearly different from the 
case of non-magnetized NSs, where there is a clear separa- 
tion in two scenarios: standard cooling (low mass) and rapid 
cooling (M > 1.6M(T) in our model) discu s sed extensively b y 
other authors (|Yakovlev fc Pethicklliooi : iPage et all 120061 ). 
In magnetized NS the fast cooling scenario is masked by 
magnetic heating, becoming hard to distinguish whether or 
not a fast neutrino emi ssion process is active in high mass 
stars, as pointed out in lAguilera et al" (2008a). This raises 
an important issue: disentangling the magnetic field initial 
distribution function is needed before we can actually con- 
strain other physical parameters that affect neutron star 
cooling. For this purpose, we use a population synthesis tech- 
nique that accounts for joint statistical properties of a set 
of objects, rather than trying to fit individual objects with 
particular cooling models. 

As in the case of the cooling scenario without extra 
heating by magnetic field decay, objects with the smallest 
masses are hotter in average. In the mass distribution func- 
tion we use they are very abundant and, since they are more 
easily detectable because of their larger thermal luminosities 
(for the same field range), most of the observed objects must 
be low mass NSs. For NSs with M > 1.3M the effect of 
changing mass is controversial, since it strongly depends on 
the assumptions about neutrino fast cooling processes. Com- 
parison between the left and right panels of Fig[T]shows that 
changes in magnetic field by a factor of ~ 3 are more im- 
portant than changing from I.3M0 to 1.76M0 for a fixed 
magnetic field. This implies that population synthesis stud- 
ies are more sensitive to varying the initial magnetic field 
distribution than varying the NS mass distribution, unless 
the population of low-field (< 10 13 G) compact objects is 
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Figure 1. Left: Comparison of cooling curves of an M = 1.25Mq NS with different values of the initial magnetic field strength. From 
bottom to top: B = 3 X 10 12 , 10 13 , 3 X 10 13 , 10 14 , 3 X 10 14 , 10 15 , and 3 X 10 15 G. Right: Comparison of cooling curves with B = 10 14 
G and different masses: M = 1.10 (top solid line), 1.25 (dots), 1.32 (dashes), 1.40 (dash-dot), 1.48 (dash-triple dot), 1.60 (long dashes), 
and 1.70 Mq (bottom solid line). 
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Figure 2. Luminosity as a function of the initial magnetic field 
strength for a M = 1.25M© NS at different ages (10 5 , 5 X 10 5 , 
and 10 6 years). 



discussed. Normally, one would expect that isolated, very 
cool objects are low field and high mass NSs, so these would 
be the best candidates to test fast cooling mechanisms with- 
out the complications due to the presence of strong magnetic 
fields. 

When we compare our model calculations with observa- 
tions, we confront emission properties of NSs (observed vs. 
calculated fluxes or luminosities), not directly their fields, 
but there is observatio nal evidence for the correlation of 
these two magnitudes (|Pons et al.ll2007l ). In the next sec- 
tion, when we establish constraints on the initial magnetic 
field distribution, the reader must keep in mind that these 
are model dependent constraints. What observations actu- 
ally constrain is the number of luminous objects, which we 
will translate to magnetic field strength using our theoretical 
models. As an example, in Fig.[2]we plot the luminosity as a 
function of the magnetic field strength for a NS model with 
M — 1.25M0 and at different ages. While the luminosity of 
young (< 10 5 yr) objects is less dependent on the magnetic 
field strength and relatively high (and therefore young NSs 



are more easily observed), the luminosity of middle aged NSs 
depends strongly of the magnetic field strength. For a given 
NS model of a certain age, the luminosity increases sharply 
above a certain value of the initial magnetic field, thus mak- 
ing more magnetized objects more likely to be observed. We 
remind again that the :r-axis indicates the initial value of 
the magnetic field at birth, not the corresponding value at 
a given age, which is always smaller as discussed above. 



3 LOG N - LOG S DISTRIBUTION FOR 
NEARBY COOLING NSS 

Previously we performed several calculations of the Log N 
- Log S distribution for NSs in the solar proximity for dif- 
ferent sets of cooling curves. Here we use the model of ther- 
mal evolution described above which includes the magnetic 
field decay. The main motivation is related to the fact that 
in the standard picture magnetic fields of magnetars decay, 
and for some of close-by NSs (the M7) there are indications 
that their fields are ~ 10 13 ' 5 G, so additional heating due 
to field decay can be important. Magnetic fields of some of 
the M7 sources are estimated by two methods: spectroscopy 
and spin-down rate. The first is based on the unconfirmed 
hypothesis that wide dep ressions in their sp ectra are due 
to proton cyclotron lines (|Haberl et al.l [2004J ) . The second 
method applies the usual magneto-dipole braking formula, 
but this is only possible when a measure of the spin period 
derivative is available. Estimates due to both methods pro- 
vide more or less consistent results within a factor of few. In 
our model with magnetic field decay, the typical strengths 
of about fewx 10 13 G are reached within fewx 10 5 yrs for ini- 
tial field s ~ 10 14 G. consis tent with age estimates for these 
NSs (see lPage et alj l|2009l ) and references therein). Temper- 
atures and spin periods of the M7 sources within our model 
are also consistent with such ages (see Sec. 6). 

To calculate the Log N - Log S distribution (number of 
objects N with a flux above S) of close-by isolated cool- 
ing NSs we use the Monte Carlo code dev eloped before 
l|Popov et alj|2003t 120051 : iPosselt et ai]|2008l ) that builds a 
synthetic population of nearby isolated NSs. A general de- 
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sc ription of the populat i on sy nthesis technique can be found 
in lPopov fc Prokhorovl (|2007). The main ingredients of the 
present model are: the spatial distribution of NS progenitors, 
the interstellar medium (ISM) density distribution needed to 
calculate the observed flux, the NS mass distribution func- 
tion, a set of cooling curves, discussed in Sec. 2, and the ini- 
tial magnetic field distribution. We now comment on each 
of this points before discussing our results. 

Progenitor s of NSs are distrib uted according to the 
model used by Posselt et al.l (2008). The contribution due 
to close-by OB-associations (the Gould Belt) is crucial. As 
before we consider the region up to 3 kpc from the Sun. 
Inside 500 pc we take as the distribution of progenitor s the 
distribution of massive stars from HIPPARCOS data (ESA 
1 19971 ). Outside this volume, most of NSs (243 per Myr out 
of the total number of 270 per Myr born inside 3 kpc) orig- 
inate in one of the OB associations. Others are distributed 
in the exponential Galactic disc. 

In our model for Log N - Log S calculations, the NS 
formation rate (in units of object per pc per year) is dif- 
ferent at different distances from the Sun due to the Gould 
Belt contribution and the non-uniform distribution of OB 
associations. Therefore, it is not straightforward to extrapo- 
late the local rat e to obtain the to t al Ga lactic rate of SNae. 
Using data from iTammann et al.l (I1994T ) we estimate that 
the value we use (~ 250 NS inside 3 kpc per Myr) can be 
rescaled to ~1.2 NSs per 100 years for the whole Galaxy. 
Here the number 250 is obtained by subtracting the addi- 
tional contribution due to the Belt (~ 2/3 of all nearby NSs 
are produced in 600 pc around the Sun). We stress that the 
value 1.2 per 100 years is only a rough estimate. Sources 
beyond ~ 1 kpc from the Sun do not contribute much to 
the Log N - Log S distribution, and in this region the NS 
formation rate in our model is twice larger than inside 3 kpc. 

As we are mostly interested in the Log N - Log S behav- 
ior at ROSAT count rates > 0.01 - 0.1 cts s" 1 , the sources 
at < 1 kpc dominate (those born in the associations forming 
the Gould Belt and in other not very far away OB associa- 
tions). That is why the global Galactic distribution (Galactic 
arms, etc.) of NS progenitors is not important for our study 
with applications to ROSAT cooling NSs. 

For the ISM distribution we used an analytical de- 
scription, which was demonstrated to be successful before 
jPopov et al.ll2003l ; |Posselt et alj200Sl ). Since here we do not 
intend to produce an accurate full-sky map of the distribu- 
tion of sources, and for computational limitations, we do not 
consider more detailed models for the ISM 3D distribution. 

For the mass distribu t ion w e use one of the variants, 
presented in IPosselt et al. |2008[) . The spectrum is derived 
using HIPPARCOS ilESAlll997l) data about clo s e-by massive 
stars and calculations bv lWooslev et all l |2002h : lHege~ etajj 
|2005l ). We use eight mass bins (1.1, 1.25, 1.32, 1.4, 1.48, 1.6, 
1.7, 1.76 M ). The first two bins contribute ~ 30% each. 
The last two less than 1%. According to this distribution 
90% of NSs are born with M < 1.45M . Such mass spec- 
trum is in agreement with mass measuremen ts of the sec - 
ondary components in double-NS binaries, fsee lStairs! (2008) 
and references therein). These NSs never accreted and they 
can be accepted as good sources for initial mass determina- 
tion (unless s o me ef fects of binary evolution are crucial). In 
IPosselt et all (|2008l ) it was shown that realistic manipula- 



tions with the mass spectrum do not influence Log N - Log 
S distributions significantly. 

We accurately integrate spatial trajectories us- 
ing the bi-modal Maxwe llian kick velocity distribution 
|Arzoumanian et al.l l2002h and the potential tr aditionally 
used in papers on isolated NSs starting with iPaczvnskil 
(1990). Nevertheless, the velocity distribution of NSs and 
the Galactic potential are not important ingredients for the 
results shown in this section, as we deal with young sources 
which do not move significantly from their birth places, and 
their velocities are nearly constant during this time. We 
tested several velocity distributions including the double- 
side expo nential with the mean ve l ocity 380 km s" 1 pro- 
posed by iFaucher-Giguere fc Kaspl <|2006h . Our results are 
not sensitive for variation between different velocity distri- 
butions which were proposed during last years and which 
successfully explain the radio pulsar observations. Also, ve- 
locity by itself does not influence the observability of sources 
under study (in contrast with, for example, studies of iso- 
lated accreting NSs). 

We calculate ROSAT counts for the Log N - Log S 
plots assuming that the local emission (there is a surface 
anisotropy consistent with the magneto-thermal evolution 
models) is purely blackbody, i.e. we neglect any non-thermal 
contribution and we do not consider effects of composition or 
magnetic fields in the atmosphere. This is partly justified by 
the fact that the M7 dominate the sample, and their non- 
thermal emission is negligible (|Haberll 120071 ) . Atmospheric 
and magnetospheric models can change the spectral energy 
distribution (not the total luminosity) significantly, resulting 
in differences in the observed flux at the detectors. However, 
we think that our sample is too small, and our knowledge 
about average properties of NS atmospheres is not mature 
enough to be included in a population synthesis scenario. 
We suspect that taking this effects into account would not 
change the results significantly, but in future works, with 
better knowledge of NSs properties and a larger sample 
it will be necessary to include these effects explicitly. We 
must also mention that we considered only heavy element 
envelopes in our model. The e f fect o f an accreted H-He en- 
velope was discussed in iPaeel (jl997h and it will be worth 
exploring in future works. 

Sources are observable in soft X-rays while they are hot. 
In our previous studies we used cooling tracks till an object 
cools down to 100 000K. In this paper we put the limit at 
300 000K due to computational reasons. However, we tested 
that this modification does not change our results, as we 
confront our results with observations of relatively bright 
sources. 

Results of the population synthesis modeling are ap- 
plied to the ROSAT all-sky survey, which is the most 
complete survey available at the energy range of interest 
(0.1-1 keV) for our research, and the most uniform sam- 
ple of close-by young objects suitable to be used in pop- 
ulation synthesis calculations. For each NS "observed" in 
the simulation we calculate ROSAT PSPC counts per sec- 
ond (S). Log N - Log S distribution for these sources is 
shown in Figs(3] [5] Filled symbols correspond to additions 
of one of the M7 NSs, empty symbols correspond to other 
sources - Vela, Geminga, B1055-52, B0656+14, and 3EG 
J1835+5918 (the second Geminga). Note, that for the last 
five sources we plot total ROSAT counts (i.e. blackbody plus 
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Figure 3. Log N - Log S distributions for different initial mag- 
netic fields. From bottom to top: 3 X 10 12 , 10 13 , 3 X 10 13 , 10 14 , 
3 X 10 14 , 10 15 , 3 X 10 15 G. "BSC" - is an upper limit from the 
ROSAT Bright Source Catalogue. Filled symbols correspond to an 
addition to the distribution one of the M7 sources, empty sym- 
bols correspo nd to an addition of PSRs (including "the secon d 
Geminga" , see iMirabal fc Halpernl ll200ll1 ; lHalpern et~aTl J2007fi ) . 

non-thermal magneto-spheric emission). If we include only 
the thermal component, then the shape of the observed Log 
N - Log S distribution is slightly changed, but clearly the 
first and the last filled points will remain as they are, and 
the change in the shape is not significant. Also, non-thermal 
contributions to luminosit ies of these PSRs are not large 
|Becker fc Truemperll 19971 ) . Error bars represent poissonian 
statistical errors (square root of the number of sources). 
"BSC" is t he upper limit fro m the ROSAT Bright Source 
Catalogue (|Voges et al.l ll999T). In the second Log N - Log 
S plot we also show as a horizontal line the most recent 
limit at 90% confidence level: TV < 31 for soft, non- variable 
NSs, and < 46 for all NSs, including hard and variable 
sources (M. Turner et al., in preparation, R. Rutledge, pri- 
vate comm.). We do not include magnetars in this sample. 
If we had considered further regions of the Galaxy (birth 
place beyond 3 kpc from the Sun), we should have included 
the additional contribution due to young magnetars, which 
can be visible from very large distances (all known young 
active magnetars are at distances > 3 kpc, according to the 
McGill SGR/AXP online catalogue 0). 

In Fig(3]we show seven curves, each one calculated for 
a single value of the magnetic field at birth (i.e. all NSs in 
the modeled population have the same field). This illustra- 
tive graph demonstrates that, for our NS model, low- field 
NSs (B < 3 x 10 13 G) cannot explain the observed sources. 
If all NSs were born with the same initial field, it should 
be in the range 3 X 10 13 — 10 14 G. Larger initial fields re- 
sult in hotter NSs, and therefore in a large number of de- 
tectable sources. Here we focus on close-by NSs with fluxes 
above ~0.1 ROSAT PSPC counts per second, as in this 
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Figure 4. \ 2 (for 10 degrees of freedom) contour levels of the 
fits to the Log N - Log S curves in the two-parameter space 
log(Bo/[G]) and °logBo' hashed line corresponds to 8% of 
NSs born with log(B /[G]) > 14.2. 

range identification of NSs among ROSAT so urces is consid- 
ered t o be mostly complete (see, for example. ISchwope et al.l 
1 19991 : ICagnoni et al.l [20021 and references therein). This is 
also confirmed by id e ntification efforts l ike th ose by e.g. 
lAgiieros etafl (120061 ); IChieregato et al.1 (|2005l ) and refer- 
ences therein. At lower fluxes (< 0.1 cts s _1 ) many sources 
might be non-identified, yet. The "second Geminga" (the 
source with the smallest count rate) is 7-ray selected, and 
this point must be taken as a lower limit for the Log N - Log 
S distribution at this flux range. Clearly, multi- wavelength 
studies (cross correlation between catalogues obtained in dif- 
ferent bands) are necessary to find more close-by cooling 
NSs using the existing data, although one of the most cru- 
cial issues for identification is probably the X-ray positional 
accura cy. Efforts in this direction, such as in lRutledge et all 
(2008), will be very useful. 

From results in Fig[3]it is clear that, unlike in previous 
works, now we have also to worry about the initial magnetic 
field distribution (B-distribution), because now our cooling 
curves depend not only on masses of NSs but also on their 
fields, and the latter effect seems to be more important. 
For most of the paper, we have chosen the birth magnetic 
field to satisfy a log-normal distribution with central value 
x c = log Bq and standard deviation cr\ og b because this type 
of distribution reproduces the observed distribution of PSRs. 
The probability of a NS to be born with a magnetic field in 
the range between Bi and B2 is then 




where x = logB and erf(i) is the error function. 

In Fig. |5]we show contour plots of the \ 2 distribution of 
fits to the Log N - Log S curves in the two-parameter space 
log Bo and <Ti og s , where Bo is given in Gauss. The best fit to 
the observational data (cross) obtained with the IDL proce- 
dure CURVEFIT is log(B /[G]) = 13.14 and cr logSo = 0.74 
but it is clear that the results of the fit are highly degener- 
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Table 1. Magnetic field distributions used in this work. Three of them (G1,G2,G3) are log-normal distributions defined by its central 
value (x c ) and width (ci gB)- Numbers in the table indicate the fraction of NSs in a given B-field bin centered in logS = x c . The other 
distributions (A1,A2, No mag) are discrete distributions and the numbers indicate the fraction of NSs with that particular field strength. 
The right column defines the linestyle used for each distribution in the plots. 
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ate, and any pair of parameters located in the central diag- 
onal band are allowed. In the figure, we also show the line 
that corresponds to the family of log-normal distribution 
functions with an 8% of NSs born with log(B /[G]) > 14.2 
(Bo > 1.6 x 10 14 G), which interestingly is nearly paral- 
lel to the band of lowest x 2 . We should stress again that 
this is a model-dependent result. For our given NS evolution 
model, this means that actually Log N - Log S curves are not 
constraining independently the average field at birth or its 
dispersion (the birth parameters determine the whole evolu- 
tion) , but simply the number of NSs born as high- luminosity 
objects, which according to our underlying physical cooling 
model is correlated to the number of magnetars. For other 
cooling models not drastically different from ours, any rea- 
sonable field distribution with approximately 8-10% of NS 
born as magnetars should also be acceptable. This is con- 
sistent with the fact the we do not observe any nearby (< 
3 kpc) magnetar. It is important to note that, since magne- 
tars are visible for a long time from large distances, results 
are so sensitive to the addition of more magnetars that even 
with ~10 observed sources we can place constraints on the 
fraction of magnetars. 

In Fig[5] we present Log N - Log S curves for six B- 
distributions (see the Table). Three log- normal distribu- 
tions (Gl, G2, G3), as they are selected from the best fit 
for Log N - Log S, pass through the observed points. For 
comparison, we also use several other variants of the B- 
distribution, summarized in the Table. Values in the Table 
correspond to fractions in each magnetic field bin normalized 
to unity. The first one is an extreme case with no NS with 
fields above 10 13 G. The other two are "hand-made" distri- 
butions (Al and A2) in both of which 1/2 of NSs belong to 
the PSR-range fields, and the rest is distributed among high- 
field objects. The curve for model A2 (dotted line) shows 
that a 30% of NS with magnetic fields > 10 14 G is already 
in contradiction with observations of the local population 
of cooling NSs. Addition of more NSs with very large initial 
magnetic fields (i.e. model Al) largely overpredicts the num- 
ber of nearby objects detected as bright thermal sources. 

In the presented graphs all Log N - Log S distributions 
are plotted for 5000 calculated tracks, each of which is used 
for all eight masses, and all data along a track is used with a 
time step 10 4 yrs. So, the results are significantly smoothed. 
We made a few additional runs for realistic numbers of NSs 
(810 NSs born during 3 Myrs in the calculated region up to 
3 kpc from the Sun). Of course, such Log N - Log S distribu- 
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Figure 5. Log N - Log S distributions for six variants of B- 
distributions. From top to bottom: Al, A2, G3, G2, Gl, No mag. 
See the Table for description of models and curve styles. In this 
plot we also added a horizontal dashed lines which correspond to 
46 and 31 sources (see the text). This is an upper limit by M. 
Turner et al. (in prep.; B. Rutledge private communication). 



tions are much more noisy. However, statistical fluctuations 
cannot change our qualitative results. For example, poisso- 
nian error bars for the data points typically bound curves 
for Gl, G2, and G3 distributions. Curves for magnetic field 
distributions with a significant amount of magnetars cannot 
explain the data even taking into account statistical fluc- 
tuations (unless, of course, a very rare strong fluctuation 
happened). 

Having in mind that observational data can be fitted 
by different field distributions, and that there are other im- 
portant parameters not explored (starting from superfiuid 
gaps and ending with properties of the ISM), the impor- 
tant message is that we can reproduce the data on the Log 
N - Log S plot with realistic distributions, a normal NS 
model, and without fine tuning. This gives us an opportu- 
nity to put a constraint on the number of magnetars. For our 
magneto-thermal evolution model, initial field distributions 
with more than 30% of NSs with initial fields above 10 14 G 
(even if all others have low fields) can be ruled out. Assum- 
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ing a log-normal initial field distribution, this also implies 
that the fraction of NSs born with B > 10 15 G can hardly be 
above a few percent. Varying other parameters of the model 
can modify this conclusion, but not dramatically, unless we 
allow for significant variations of the birth rate. 



4 LOG N - LOG L DISTRIBUTION FOR 
GALACTIC MAGNETARS 

Motiv ated by a recent important paper by iMuno et al.l 
(2008), and as a consistency check of the results of the previ- 
ous section, we also consider a simple model for the Log N - 
Log L distribution of highly magnetized NSs. In lMuno et al.l 
(2008) the authors analyze a large set of Chandra and XMM- 
Newton data (nearly one thousand exposures) to search for 
new magnetars looking for pulsating sources in the range 5- 
200 seconds. No new candidate was found, but this fact could 
be used to place upper limits on the total number of Galactic 
magnetars with different properties. These authors estimate 
that the number of magnetars with L > 3 x 10 33 erg s _1 and 
pulsed fraction larger than 15% is below 540, and the num- 
ber of easily-detectable magnetars is 59lj2 . Of course, some 
(perhaps many) magnetars can have lower pulsed fractions 
and not bei ng detected as X-r ay PSRs. In this respect the 
numbers bv lMuno et all (2008) are representative of a frac- 
tion of the total population, and total limits maybe larger 
than claimed. These are, to our knowledge, the best obser- 
vational limits on the number of such objects. On the other 
hand, the fact that we observe some magnetars can be used 
to place lower limits. These limits are too weak to favor 
one model against another, but it is illustrative to compare 
some of the models that best fit the close-by, intermediate 
field NSs to the magnetar population. 

In our model the luminosity is of thermal origin, because 
we do not consider magnetospheric processes that may cause 
the non-thermal emission. However, since in magnetars both 
thermal and non-thermal components are supposed to be 
powered by magnetic field decay, and energy conservation is 
satisfied by construction in our evolutionary models, we can 
conclude that the total energy release is calculated correctly, 
although the spectrum can be different. 

We must stress that in this part no Monte Carlo sim- 
ulation is done. Instead, we use complete cooling tracks of 
NSs with different masses and magnetic fields to estimate 
the whole Galactic population of NSs with a given luminos- 
ity. Absolute numbers are obtained by normalization to the 
total birth rate of NSs. We use the Galactic NS formation 
rate equal to 1/30 yr s -1 . This is close to the upper limit 
for NS formation rate lKeane fc Kramer! l|2008l ). The uncer- 
tainty in the NS format ion rate is a factor ~ 2 — 3 (see also 
iKeane fc Kramerl ()2008h and references therein) which may 
shift curves in Fig(6] This NS formation rate is not directly 
related to the rate used for the Log N - Log S calculations 
above. In the case of close-by cooling NSs the rate of NS 
formation is determined by the properties of the Gould Belt 
and close- by (< 3 kpc) OB-associations. Here, in the case 
of magnetars, we are interested only in the gl obal galactic 
rate o f NS formation (see a similar discussion in lGill fc Hevll 
120071 '). 

As we calculate distribution in luminosity, not in flux, 
we do not take into account interstellar absorption. Calcula- 



tions are done for the same cooling curves, B-distributions, 
and mass distribution as used for Log N - Log S calcu- 
lations. Formally, very nearby sources with low fields can 
also contribute to the observed Log N - Log L distribution, 
but in the range of luminosities we are interested in (> few 
xlO 33 erg s _1 ) their contribution is negligible. 

In Fig(6] we show Log N - Log L distributions calcu- 
lated for the same six B-distributions describe d in the Ta- 
ble. W e compare our curves with the data by IMuno et al.l 
(2008): the upper limit to the number of magnetars with 
L > 3 x 10 33 erg s -1 and pulsed fraction larger than 15% is 
540 and the number of easily-detectable magnetars is 59^t§2 • 
Since Muno et al. (2008) consider as easily detectable mag- 
netars two types of objects: bright magnetars with small 
pulsed fraction, and dim, but with very large pulsed frac- 
tion, we show this as a rectangular region through which 
satisfactory models should pass. As for real observations, 
at the moment we know 5 SGRs (plus candidates), and 
10 AXPs (plus transient objects and candidates, see the 
McGill group on-line catalogue). These are also shown in 
the figure (diamonds) for comparison. In the on-line cata- 
logue the luminosity is given for the range 2-10 keV. Some 
magnetars also de monstrate significant hard X-ray emission 
(Mereghctti 2008). It is not included in the plot, but in log- 
scale the shift is not crucial. Note that this sample is not 
complete, so it must be taken as a lower limit for the pre- 
diction of theoretical models. At the very bright tail it is 
possible that the sample is close to complete, because there 
should be no more very bright L > 10 35 erg s -1 magnetars. 

Modeled Log N - Log L distributions start to flatten 
at L ~ 10 33 erg s _1 . This value corresponds to the weak- 
est known magnetars. The results from the previous section 
show that all acceptable log-normal distributions for middle- 
aged thermally emitting NSs are similar in this luminosity 
range. They all predict about a thousand NSs above this lu- 
minosity. Note that we assume that all sources are persistent 
and the models considered in this work only include steady 
mag netic field Ohmic decay. In addition to purely Ohmic de- 
cay, iPons fc Geppert] (|2007h found that the Hall drift may 
contribute noticeably to accelerating the dissipation of mag- 
netic fields in young NSs. This Hall phase lasts a few 10 3 -10 4 
years and is characterized by an intense exchange of mag- 
netic energy between the poloidal and toroidal components 
of the field and by the redistribution of magnetic field en- 
ergy between different scales. It can be expected that such 
rearrangements and the relatively rapid field decay may en- 
hance the average luminosity and result in crustal breaking 
and active stages (bursts, flares), as can be observed in mag- 
netars. This would increase the number of bright sources 
with L ~ 10 35 erg s _1 , and this expectation is preliminary 
confirmed by some artificial models in which we tried to 
take into account the possibility of this transient behavior 
in young highly magnetized NSs. Alternatively, one can also 
consider a fraction of NSs with larger internal toroidal fields 
that have larger luminosities, but this introduces yet another 
free parameter in the problem. A careful study of the first 
few thousand years of an ultra-luminous magnetar's life and 
its transient epochs is out of the scope of this paper, that 
focuses on long-term evolution and statistics. 
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Figure 6. Theoretical Log N - Log L curves for Galactic magne- 
tars compared with observational data and constraints. Diamonds 
show the Log N - Log L distribution for known magnetars (from 
the McGill group on-line catalogue), the square indicates the limit 
of 540 weak AXP (Muno et al. 2008), and the box corresponds to 
the estimated number of "easily detectable magnetars" (Muno et 
al. 2008). We use the same linestyle for B-distribution models as 
in Fig|5](see the Table). 



5 EVOLUTION OF PSRS AND THE P - P 
DIAGRAM 

We turn now to PSRs. We have performed Monte Carlo 
simulations to generate a synthetic PSR population and 
confront our models with observations. The methodology 
employed in the simulati o ns cl osely follows the work by 
IFaucher-Giguere fc KaspH (I2006T I but some parameters of 
their model are allowed to change according to the results 
of our previous sections. The main goal of this section is 
to answer the following question: can we obtain a syn- 
thetic PSR population compatible with the observed one 
and consistent with our previous description for magnetars 
and close-by isolated neutron stars ? To this end, we start 
from the optimal populati o n mo del parameters obtained by 
IFaucher-Giguere fc K aspi ( 200(3) and modify only the ini- 
tial period and magnetic field distributions to account for 
the effect of magnetic field decay consistent with our model. 

To generate the PSR synthetic population we first 
choose the parameters o f the NS at birth closely fo l lowin g 
the model described by IFaucher-Giguere fc KaspH |2006), 
which we briefly summarize. The age of the NS is chosen 
randomly in the interval [0, t ma x], where t ma x = 500 Myr. 
This is shorter than the age of the Galactic disk, but it is 
enough for our purpose because PSRs older than this age 
have crossed the death line (assuming standard magnetic 
dipole braking) and are no longer visible as PSRs. The place 
of birth is obtained according to the distribution of their 
progenitors (massive Population I stars) which are mainly 
populating the Galactic disk and more precisely its arms. 



The velocity at birth is distributed according the exponen- 
tial distribution with a mean value of 380 km s _1 .[j 

The spin period of the star at birth, Po, is chosen from 
a normal distribution with a mean value of < Po > and 
standard deviation ap . Of course, only positive values are 
allowed. The initial magnetic field at the magnetic pole is 
obtained from a log-normal distribution with mean value 
< log(£>o/[G]) > and standard deviation o\ os b ■ 

Once we have chosen the properties of the NS at birth 
we solve the appropriate differential equations to obtain the 
position, period and magnetic field at the present time. We 
use a smooth model for the Galactic gravitational pot ential 
(|Kuiiken fc Gilmord fl989l ; ICarlberg fc Innanenlll987l ). The 
period evolution is obtained by assuming that the rotation 
energy losses are due to magnetic dipolar emission (orthogo- 
nal rotators), where the magnetic field is obtained from our 
magneto-thermal evolutionary models described in section 
2. 

At the end of the Monte Carlo simulation we end up 
with a synthetic population of PSRs to be compared with 
a given observed sample. We use the PSR s detected i n the 
Parkes Multibeam Survey (PMBS) sample l|Lvnell2008l ) and, 
to limit the contamination of our sample by recycled PSRs, 
we further ignore the PSRs with P < 30 ms or P < 0, and 
those in binary systems. With this restrictions, our result- 
ing sample contains 977 objects. We use the parameters for 
detectab ility in the survey, radio l umino sity and beaming 
given m IFaucher-Giguere fc Kaspl (|2006l ). For comparison, 
The P — P diagram for the sample of real pulsars retained 
in our analysis is shown in Fig. [7] together with some evolu- 
tionary tracks of the model used in the analysis. For lower 
fields, NSs move nearly along constant magnetic field lines. 
We must remark again that, in magnetars, a somewhat en- 
hanced field decay is expected to happen during the initial 
Hall stage. This non-linear term is not yet included in our 
numerical simulations and we expect a more vertical initial 
trajectory for objects in the upper right corner of the dia- 
gram. 

The problem of defin ing an optimal model has b e en we ll 
discussed in Sect. 3.7 of Fau cher- Giguere fc Kaspil l|2006(h 
and we have adopted their approximate approach that re- 
quires of some human judgment, rather than attempting to 
cover a huge parameter space with our limited computa- 
tional resources. We have explored the parameter space in 
the region that best fits the properties of thermally emit- 
ting NSs and magnetars, as described in previous sections. 
Without pretending to perform a rigorous fully quantita- 
tive analysis, we have considered the Kolmogorov-Smirnov 
goodness-of-fit test to quantify our statistical analysis. 

In Fig. [8] we show P — P diagrams for typical Monte 
Carlo realizations and their corresponding distributions of 
observed PSR periods and magnetic fields (averaged over 50 

3 From the point of view of velocity and initial spatial distribu- 
tions these assumptions differ from those used for Log N - Log S 
calculations. However, as we describe in Sec. 3, Log N - Log S cal- 
culations are not very sensitive to the velocity distribution and to 
large scale spatial distribution. The first is due to relatively small 
ages of studied sources. The second, because in Log N - Log S 
at significant ROSAT count rates close-by sources dominate. The 
same can be said about differences in the models for the Galactic 
potential in P —P and Log N - Log S calculations. 
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Figure 7. P — P diagram for the sample of 977 real pulsars re- 
tained in our analysis (small filled circles) with seven evolutionary 
tracks with different initial magnetic fields. For comparison, we 
also show AXPs and SGRs (open diamonds) three of the M7 (tri- 
angles), and RRATs (filled diamonds). For RRATs new data from 
McLaughlin et al. (2009) is used. Color of the tracks reflects the 
NS temperature, and long dashed lines indicate true ages of 10 4 , 
10 5 and 10 s yr. The color bar shows temperature in logarithmic 
scale. 



realizations of 977 detectable PS Rs). The upper and lower 
panels show the optimal models of lFaucher-Giguere fc Kaspil 
( 2006) and this work, respectively. The values for the distri- 
bution of the optimal model assuming there is no field decay, 
are < log(S /[G]) >= 12.95 and cri ogBn =0.55, < Pp >= 0.3 
s, and op Q — 0.15 s. Note that IFaucher-Giguere fc Kaspil 
(2006) use the magnetic field at the equator and the value 
they give is < log(B /[G]) >= 12.65. 

Among all the models analyzed, we find that our op- 
timal model with realistic magneto-thermal evolution of 
NSs corresponds to < log(i?o/G) >= 13.25, o"i ogB() =0.6, < 
Po >= 0.25 s, and ap a = 0.1 s. Since our model includes field 
decay, we find that our average ini tial magnetic field is about 
a fact or of 2 larger than that of IFaucher-Giguere fc Kaspil 
(2006), and the distribution is also slightly wider. Because 
of this larger average field, the initial period distribution 
has to be shifted to lower values to obtain a visible syn- 
thetic PSR population statistically similar to the observed 
population of PSRs. 

Our main conclusion is that, within the parameter space 
that best fits the observed population of nearby, thermally 
emitting NSs and magnetars, we can also find an optimal 
parameterization that satisfactorily explains the observed 
PSR population. Therefore, it is possible to describe simul- 
taneously different families of the neutron star zoo with a 
single underlying physical model. Interestingly, a combined 
statistical analysis of PSRs and thermally emitting NSs al- 
lows to break the degeneracy in the parameter space that 
arises when we try to work with a single family. Models with 
< log(B /[G]) >= 13.3 or < Iog(B /[G]) >= 13.2, or with 
wider or narrower distributions (cri ogBo =0.7 or o"i ogSo =0.5) 
give a much worse result in the Kolmogorov-Smirnov test, 



obtaining in all these cases P-values < 0.01 s. Therefore, 
having fixed a NS model and initial magnetic field geome- 
try, and within the confidence region obtained in Log N - 
Log S analysis, we conclude that the observed distribution 
of radio PSRs is only consistent with values in a narrow 
vicinity of the optimal model. In the future, as we improve 
the NS evolutionary models results may change, but the in- 
teresting fact is that a combined analysis turns out to be 
very restrictive and breaks the degeneracy obtained in the 
study of populations of only nearby thermally emitting NSs 
or Galactic magnetars. We leave a more extensive study of 
the influence of these or other parameters (velocity distri- 
bution, birth rates, etc.) for future work. 



6 DISCUSSION AND FINAL REMARKS. 

In this paper we presented a multi-component population 
synthesis study. The final goal in this approach would be 
to make a complete population synthesis with a unique NS 
physical model that consistently explains all known types 
of young NSs just varying parameters such as the NS mass, 
age, or the strength and geometry of the magnetic field. Still, 
even after more sophisticated theoretical calculations are 
available, comparison with the data will proceed by steps, 
confronting each piece of simulated data with observational 
data of some population, in a similar way to this paper. 
One reason for engaging multi-population studies can be il- 
lustrated as follows. From Fig. 3 it is visible that the Log N 
- Log S distribution can be explained by a single field model 
(with moderate additional heating, in our framework). How- 
ever, the Log N - Log L distribution for magnetars, as well 
as the P-P plot for radio pulsars, of course, cannot be ex- 
plained by this model for different reasons. In the first case 
case, because no bright magnetars will be observed; in the 
latter, because it will be impossible to reproduce the main 
part of the pulsar population. 

To summarize, we believe that the approach we use is 
well motivated by a necessity to have a natural model with- 
out ad hoc assumptions about different subpopulations. Note 
that here we tested a very "smooth" model, in a sense that 
we do not put by hand initially distinct populations (magne- 
tars, M7, PSRs, etc.). In each part of our study (Log N - Log 
S for the M7, Log N - Log L for magnetars, P-P for PSRs) 
we model "just NSs" using the same initial magnetic field 
distribution (which in our model is the main parameter), 
without specifying unique particular feat ures for a given sub- 
population, as it i s usually done (see e.g.lPopov et alj 120061 ; 
iGill fc Hevlll2007l ; iKeane fc Kramerll2008l ). 

For example, XDINSs are not a separately defined class 
with initially distinct properties, but they just appear as a 
population coming out from a smooth unique initial distri- 
bution, and with the following features: 

• The magnetic field in these sources has significantly de- 
cayed from larger initial values, so no magnetar-like activity 
is present; 

• Due to their large initial fields spin periods are long, so 
no radio pulsar activity is observed (may be due to narrow 
beams) ; 

• They are still young enough for th e magnetic field t o 
be still decaying (< 10 6 yrs, according to lPons et alj (|2009h ) 
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Figure 8. P — P diagram for typical Monte Carlo realizations and distributions of observed PSR periods and magnetic fields. The 
distributions show the average of 50 realizations (error bars indicate standard deviatio n) compared to the PMBS sam ple distribution 
(red lines). The upper panels show the results for the optimal model without field decay ( Fauchcr-Gigucrc &; Kaspi 2006), and the lower 
panels correspond to the optimal model consistent with our simulations of magneto-thermal evolution of NSs with field decay. On each 
histogram the associated Kolmogorov-Smirnov P-value is displayed in the upper right corner. 



and heating the crust, so that they can be observed as rela- 
tively bright thermal sources. 

In our model a typical M7-like source has an initial field 
of B ~ 10 14 G. The period is P ~ 7 sec at a true age of 
~ 5 10 s yrs (spin-down age ~ (8 — 9) 10 5 yrs) . At this time 
such an object has B ~ 4 10 13 G and T ~ (6-7) 10 5 K. The 
additional energy being input by field decay in these sources 
can be estimated in each case by multiplying the average 
magnetic energy density (~ 10 26 erg cm~ 3 ), by the crust 
volume (~ 10 18 cm 3 ) and dividing by the typical Ohmic 
decay timescale under such conditions (crust density and 
temperature), which is ~ 10 13 s. This gives roughly ifj 30-31 
erg s _1 available from magnetic field decay. A fraction of 
this energy can be radiated in the form of neutrinos, but 
the energy reservoir to keep the star warm is still very im- 
portant. 

We tested average ages and typical parameters (periods, 
magnetic fields, etc.) of NSs which contribute to the range 
of Log N - Log S between 0.1 and 10 cts s~ , where all 
observed sources are located. On average modeled NSs with 
inital fields 10 14 -3 10 14 G which contribute a lot to this 
range have ages 2 10 5 - 5 10 5 yrs and fields (at the moment 
of observation) ~ 7 10 13 G. Periods are distributed between 
~2-20 sec, in correspondence with observed properties of 
the M7. Of course, some fraction of lower field NSs are also 
found in this range, in correspondence with observations of 
cooling radio pulsars (Vela, Geminga etc.). Detailed study 
of the modeled population in the range 0.1-10 cts s _1 is in 
progress, and will be presented elsewhere. 

In any population synthesis approach inevitably it is 
necessary to make simplifications and to neglect some de- 
tails. Partly, simplifications made in this study are justified 
by low statistics of known sources, partly by uncertain prop- 
erties of objects under study. One of the main problems in a 
population synthesis study is related to possible correlations 
between different parameters. Some correlations are irrele- 
vant for our purposes (i.e., the correlation between direction 
of velocity and spin axis) but others can be crucial. For ex- 
ample, in our study we assume that masses, initial magnetic 
fields and spin periods are not correlated. However, this is 
very uncertain for all subclasses of NSs. If these three pa- 
rameters are correlated our results may change. For magne- 



tars it was propose d that they can have massive progenitors 
l|Muno et alj|200q ). and normally mor e massive progenitor s 
are expected to produce massive NSs i|Wooslev et al.ll2002r i. 
Thus, we can expect a correlation between field and mass for 
magnetars but not for other NS classes. Such correlations, 
valid only for a not well identified part of a population, are 
very hard to confirm or to rule out. In Fig. [1] one can see 
that an increase in the value of the magnetic field is much 
more influential (on the thermal evolution) than changes in 
mass. That is why we think that moderate mass- field corre- 
lation in the case of magnetars has very little influence on 
our results. Other effects like rotation and mass-loss would 
complicate the situation even more. 

We also neglect all possible effects related to the fact 
that a significant fraction of NSs are born in binary sys- 
tems. Potentially, this can be important for magnetars if 
their large magnetic fields are generated by amplification 
due to rapid rotation of the protoNS. The precise mecha- 
nis m of magnetar formation is still unknown. In some mod - 
els ||Popov fc Prokhorovl 120061 ; iBogomazov fc Popovl I2009T ) 
magnetars are born only in binary systems where a progen- 
itor core was spun-up due to accretion or tidal interaction, 
similar to the main scenario for gamma-ray burst progeni- 
tors. 

In our simplified standard scenario (it is clearly visible 
in FigO the M7 sources do not seem to be descendants of 
extreme magnetars, but the evolutionary link with some of 
the AXPs is possible. However, due to numerical limitations 
and the very different timescales, our model does not in- 
clude the possibility of enhanced magnetic field dissipation 
during the fast initial Hall stage, nor transient phenomena 
that change a quiet X-ray emitter into an active bursting 
source, and back. Understanding the short term violent be- 
havior of young magnetars may help to reconcile even the 
highest field objects with the M7. 

One possibility to test this hypothesis is to look at the 
velocity distribution of these types of sources. Indeed, veloc- 
ity is a good invariant on the time scale of ~1 Myr. If AXPs 
are rapidly moving objects, as it was popular to assume some 
years ago, but M7 have relatively small velocities, then one 
would conclude that the two populations are not related. 
Interestingly, there is a recent measure of the transverse ve - 
locity of the magnetar XTE J1810-197 l|Helfand et al.ll2007M . 
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The measured velocity is slightly below the average for nor- 
mal young neutron stars, indicating that the mechanism of 
magnetar birth need not lead to h i gh NS velocities. On the 
other hand, recently iMotch et all $2009) reported new ve- 
locity measurements for the M7 sources, and one of them 
(RX J1308. 6+2127) appears to be a fast object. In the near 
future we expect that the measurement of proper motions of 
some of these sources with new observations will contribute 
to our understanding of their evolutionary link, if any. 

In our study we also assumed that the magnetic field 
structure is similar for all NSs and we only rescale the nor- 
malization value for a fixed geometry. This is obviously an 
arbitrary choice. There are several issues related with the 
magnetic field geometry that need further investigation. If 
the magnetic field is not supported by currents located in 
the crust, but instead by superconducting currents in the 
core, the field will dissipate on much longer timescales, and 
the heating mechanism we assume is less important. The 
amount of energy stored in an internal toroidal field, com- 
pared to the corresponding energy of the measured dipolar 
component is also unknown, and it is an important param- 
eter that can significantly alter the results. We have chosen 
toroidal fields which maximum amplitude is a factor of 2 the 
dipole estimate (however, the energy stored on the toroidal 
field is only about a 10% because it is confined to a small vol- 
ume), in agree ment with recent studies on MHD equilib rium 
configurations (|Ciolfi et al . 20091; lLander fc Jo nes 2009), but 
this remains an open question. Some models we tried with 
toroidal fields one order of magnitude larger produced much 
hotter objects, and overpredicted the observed number of 
isolated NSs and magnetars, unless we reduce significantly 
the birth rate. 

In this paper we have not attempted to adjust our cool- 
ing curves for low initial magnetic fields in such a way that 
the local population of normal PSRs with detected thermal 
emission is perfectly explained. Future more extensive cal- 
culations exploring other parameters (superfluid gaps, crust 
physical properties, etc.) can help to fit better the popula- 
tion of thermally emitting local NSs, and a better knowledge 
of the mass spectrum can be important for low-field stars, 
too. For now, we can work with the known 7 XDINS (for 
which field decay is probably important) and 4-5 normal 
close-by PSRs (i.e. lower field sources without additional 
heating). In fact, the lowest curve in Fig[5] predicts only 1-2 
objects with flux larger than 0.1 cts s" 1 , while we know four 
near-by PSRs (Vela, Geminga, B1055-52, and B0656+14). 
The low field cooling curves are more sensitive to details of 
the equation of state, superfluid gaps, and neutrino emissiv- 
ities (and the NS mass) than high field models, and changes 
of this parameters (i.e. the neutron gap in the core) can 
shift up or down the curves. Ideally, the low field population 
should be used in a separate analysis to constrain parame- 
ters of the interior physics but, given the problem of very 
low statistics we are facing, it seems meaningless to split our 
populations in more subgroups at this stage. We made some 
tests to see if it can be done but, with the low statistics we 
manage and with the present day uncertainties about de- 
tails of NS cooling, we do not think that an extensive inves- 
tigation of microphysical parameters can contribute much 
to our understanding of NS properties. We certainly need 
more data before making strong cases in favor of particular 
cooling models. 



Despite many a ttempts (see lAgiieros et al.l 120061 ; 
IChieregato et al.ll2005l and references therein) the number of 
M7-like sources is not increasing significantly. More M7-like 
NSs can be found using deep XMM-Newton and Chandra 
observation s . A g ood example is the new source found by 
iPires et al.l l|2009h . Another possibility is to have a sample 
of 7-ray sources selected by Fermi/GLAST, but they should 
be not M7-like, but PSRs with radio beams not pointing 
towards us (similar to the "second Geminga"). The most 
promising way to increase the number of M7-like sources is 
related to the future eROSITA instrument aboard Russian 
satellite Spectrum-X-Gamma. New sources are expected to 
be dimmer, youn ger, hotter and furth er away than the seven 
ROSAT sources ijposselt et £01120081 ). and eROSITA (with 
non-truncated sensitivity at low energies) will be a per- 
fect tool for them. For magnetars there is some hope to 
have many more candidates due to the MAXI instrument 
aboa rd the International Space Station (Nakagaw a et all 
2009). Then, with increased statistics, it will be useful to 
come back to study a separate Log N - Log S distribution 
of close-by cooling PSRs. 

To summarize, the methodology we employ in this arti- 
cle looks very promising to constrain NS properties as more 
data coming from future missions arrive and more PSRs 
are found by radio surveys. A future increase of the ob- 
servational data will allow to perform much more detailed 
population synthesis and to establish evolutionary links (if 
they exist) between different classes of isolated NSs, and to 
understand better the general evolution of these sources. 
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